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ABSTRACT 

We present a hydrodynamical code for cosmological simulations which uses the Piece- 
wise Parabolic Method (PPM) to follow the dynamics of gas component and an N- 
body Particle-Mesh algorithm for the evolution of collisionless component. The gravi¬ 
tational interaction between the two components is regulated by the Poisson equation 
which is solved by a standard FFT procedure. In order to simulate cosmological flows 
we have introduced several modifications to the original PPM scheme which we de¬ 
scribe in detail. Various tests of the code are presented including adiabatic expansion, 
single and multiple pancake formation and three-dimensional cosmological simulations 
with initial conditions based on the cold dark matter scenario. 

Key words: Hydrodynamics - Methods: numerical - Large-scale structure of the 
Universe 


1 INTRODUCTION 

Most of the present cosmological models are based on the assumption that in the universe two different kinds of matter are 
present: the baryonic matter, which is directly observed and forms all of the bright objects, from stars to the hot gas present 
in X-ray clusters, and a dark, collisionless component which accounts for most of the gravitational mass in the Universe. 
The evolution of this system can only be described by treating both components at the same time, looking at all of their 
internal processes and considering their mutual interaction. In this way, making a suitable choice of the initial conditions, 
one can describe the evolution of cosmological structures and estimate all of the related physical observables. This can only 
be achieved by using numerical simulations which allow a general description of the non linear evolution of the structures. 
In particular N-body techniques (Hockney & Eastwood 1981; Efstathiou et al. 1985; Barnes & Hut 1986) have proved to be 
particularly effective for cosmological problems in which the dynamics is controlled by gravitational forces as in the case of 
the dark matter and also for baryonic structures on very large scales (more than 50/U 1 Mpc). With these methods matter is 
described as a set of collisionless particles and its dynamics is governed by the Boltzmann equation. 

The N body approach, however, is not in general suitable for describing the behaviour of the baryonic component which 
is also influenced by pressure forces, heating and cooling processes. The inclusion of all these phenomena requires an enormous 
amount of computational resources as they act on a very wide range of scales. The thermal input from gravitationally induced 
shocks is, for example, maximally effective at about 10 Mpc, while cooling processes are most important on scales less than 0.1 
Mpc. Therefore to take into account all of these processes requires a minimum of 10 6 resolution elements in three dimensions. 
The lack of adequate computational resources has delayed the development of hydrodynamic cosmological codes and only in 
recent years have a number of numerical schemes been proposed for following the evolution of baryonic matter. 

A first family of these techniques derives directly from the N-body methods. It is called “Smoothed Particle Hydrody¬ 
namics” (SPH) and represents fundamental fluid elements in terms of particles. The SPH methods are intrinsically Lagrangian 
and, following the fluid elements in their motion, have high spatial resolution and give an accurate description of high-density 
regions, where particles tend to concentrate. On the other hand they cannot treat properly low-density regions where few 
particles are present and mass resolution is poor. For a detailed introduction to these methods we refer to Hernquist & Katz 
(1989), Evrard (1990) and Steinmetz & Muller (1993). 
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Recently Gnedin (1995; see also Gnedin & Bertschinger 1996) has presented a new approach to cosmological hydrody¬ 
namics called SLH (softened Lagrangian hydrodynamics) which utilizes a high resolution Lagrangian hydrodynamics code 
combined with a low resolution Eulerian solver. Its properties are intermediate between the Eulerian and SPH approaches. 

Multi-dimensional hydrodynamic codes not based on SPH, usually adopt an Eulerian approach and dynamical equations 
are solved on a fixed (or adaptive) grid. Mean values of the fluid quantities are computed in each grid cell by solving the 
equations of conservation of matter, momentum and energy density once the equation of state for the matter is given. Eulerian 
methods have good mass resolution, and can describe low-density regions better than SPH but they are spatially limited by 
the cell size. An Eulerian approach has been adopted in the majority of the hydrodynamical codes developed for studying 
large scale structures (Chiang, Ryu & Vishniac 1989; Cen 1992; Ryu et al. 1993; Bryan et al. 1995; Quilis, Ibanez & Saez 
1996; Sornborger et al. 1996). For a comprehensive comparison between the Lagrangian and Eulerian approaches we refer to 
Kang et al. (1994). 

In this paper we describe a cosmological hydrodynamic code that we have developed using the Piecewise Parabolic Method 
(PPM) introduced by Colella & Woodward (1984). This is a higher order extension of Godunov’s shock capturing method 
(Godunov 1959; 1961). It is at least second-order accurate in space (up to the fourth-order, in the case of smooth flows and 
small timesteps) and second-order accurate in time. The high accuracy of this method allows minimization of errors due to 
the finite size of the cells of the grid and leads to a spatial resolution close to the nominal one. In a cosmological framework, 
the basic PPM technique has to be modified to include the gravitational interaction and the expansion of the universe. In 
the gravitational collapse of cosmological structures, extremely supersonic flows are generated and so the thermal energy 
cannot be computed accurately from the total energy, as is usually done. Thermal energy is only a small fraction of the total 
energy and so negative values or a spurious heating can both be found as a result of numerical errors. This difficulty has been 
overcome by computing simultaneously the total and the internal energy. 

The PPM algorithm has already been used for building a cosmological code by Bryan et al. (1995), however our approach 
differs in several respects from theirs. In the Bryan et al. work, the one-dimensional time integration is done by first performing 
a Lagrangian step and then remapping the results onto an Eulerian grid expressed in the usual coordinates comoving with 
the mean Hubble flow. We instead adopt a single-step Eulerian scheme. The construction of the effective left and right states 
for the Riemann problem is then more complicated than in the Lagrangian case, since the number of characteristics reaching 
the edge of a zone is not constant. On the other hand, this choice allows us to include in the characteristic equations both 
the gravitational interaction and the expansion of the universe and then the effect of all source terms is accounted for to 
second-order. Moreover, an Eulerian approach seems to be preferable in the case of complicated three-dimensional structures 
(Colella & Woodward 1984). 

The hydrodynamical part has been coupled to a Particle Mesh (PM) N-body code that describes the evolution of the 
dark component. The standard PM code has been modified to work with a non-constant timestep equal to that used in the 
hydrodynamical integration. The coupling is obtained by calculating the gravitational field due to both the components by 
the usual FFT procedure. 

The plan of the paper is as follows. In Section 2 we present the basic equations written in comoving coordinates. In 
Section 3 we describe the basic numerical method and the modifications needed for cosmological calculations. The numerical 
tests and results from cosmological simulations are presented in Section 4, while conclusions are drawn in Section 5. A more 
detailed description of the numerical scheme and of the method used for computing the fluxes are presented in the appendices. 


2 DYNAMICAL EQUATIONS 


The cosmological hydrodynamic equations for a self-gravitating gas in coordinates comoving with the expanding universe are: 

(1) 


do Id, , 

TT7 4 -77- (QVj ) — 0 

at a oxi 


dgvi Id . a 1 d(p 

~dt ^ alkcj {0ViVj +P5i ’ j) = —a QVi ~ a Q dP t ’ 

dE 1 <9 r , a 1 d(f> 

+ — ——[(E + p)vj] = - 2 -E - gvj- 


dt a dxn 


dxi 


( 2 ) 

( 3 ) 


where Xj and Vj are the j'-th component of the comoving coordinate and of the proper peculiar velocity, respectively; t is 
the cosmic time; a is the cosmological expansion factor; (p is the peculiar gravitational potential; q is the comoving matter 
density; E is the total energy per unit comoving volume; p is the comoving pressure. 

We treat the baryonic matter as an ideal gas with adiabatic index 7 = 5/3 and its equation of state is: 


P = (7- 1) (e- \qv 2 ^ . (4) 

where v 2 is the square norm of the velocity vector. 
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The dynamics of collisionless matter is described in comoving coordinates by the following set of equations: 

dxj 1 
~dt ~ a Vj ’ 


dvj a 1 dd> 

—A H —Vj =— 
at a a oxn 


The peculiar gravitational potential is computed by solving the Poisson equation: 
_ 47rG 

^ a \@tot Qo) 5 


( 5 ) 

( 6 ) 

( 7 ) 


where g tot is the total comoving density (baryonic plus dark matter) and g 0 is the comoving mean density. 

In the code we have normalized the variables using as basic units the final time to, the comoving cell size xo and the 
comoving mean density g 0 at time to - The expansion factor at the final time is set as ao = 1. The other variables are normalized 
as follows: Vj/(x 0 /t 0 ), E/(g 0 xl/tl), p/{g 0 xl/tl) and <t>/{xl/tl). 

The physical temperature T is related to the previous normalized quantities by 

TT w 


T = 


K b 


'Xo 
k to ) 


where Kg is Boltzmann’s constant and m v is the proton mass. 


3 CODE STRUCTURE 
3.1 The PPM Scheme 

The choice of the PPM scheme for the integration of liydrodynamical equations has been suggested by the particular problems 
which one has to face in a cosmological computation. The presence of highly supersonic flows and the formation of strong 
shocks during the collapse of structures require a very accurate treatment. Moreover the description of phenomena which act 
simultaneously on a very wide range of scales needs as much resolution as possible with a given number of cells. 

The PPM algorithm, developed by Colella & Woodward (1984), belongs to a class of schemes in which an accurate 
representation of flow discontinuities is made possible by building into the numerical method the calculation of the propagation 
and interaction of non-linear waves. There are two possible formulations of the PPM scheme: liydrodynamical equations can 
be solved either in a single Eulerian step or with a Lagrangian step followed by a remapping onto the Eulerian grid. We 
have adopted the first approach as it appears more suitable for problems with complicated spatially dependent source terms. 
A detailed presentation of the finite difference equations used in the code are presented in Appendices A1 and A2: here we 
describe the basic operations characterizing the PPM scheme. 

In the PPM method we start from the knowledge of a set of zone-averaged values of liydrodynamical quantities. To 
update these averages we solve liydrodynamical conservation equations: this requires the estimate of the fluxes of the conserved 
quantities at the zone interfaces which is done according to the following one-dimensional procedure. First, we construct a 
piecewise parabolic (third-order) one-dimensional interpolation function for each of the liydrodynamical variables and for 
the source terms. This leads to a more accurate representation of smooth spatial gradients as well as a sharper description of 
discontinuities. The interpolation polynomial w is computed such that: 

3 + 1/2 

> (9) 

- 1/2 

where Wj is the mean value of the interpolation polynomial over the j-th zone, Xj ± 1 / 2 are the positions of the zone interfaces 
and Aa; = Xj + \/ 2 — Xj_\/ 2 . It is also required that inside each zone w(£,) is continuous and monotonic. 

The interpolated distributions can then be used to estimate properly the average values of the variables to the left 
and right of each zone edge. This is performed in two steps. In the first step, one determines the characteristic domain of 
dependence with respect to a zone edge, that is the regions of space to the left and to the right of the edge that contain 
all of the information that can reach the zone edge during the current timestep. Then, using the interpolated distributions, 
the mean values of hydrodynamic variables are computed over the domains of dependence on each side of the edge. In the 
second step, we solve the Riemann problem using as initial left and right states the mean values of liydrodynamical variables 
computed previously on each side of the zone edge. The solution of this gives the time averaged values of the variables at the 
zone edge which are then used to compute inter-cell fluxes and hence the fluid variables at the new time. This closes the basic 
one-dimensional PPM loop. 

The extension to more dimensions is achieved using a direction splitting procedure (Strang 1968). If we indicate with Lk 
the integration in the fc-th direction, the updated value for the generic liydrodynamical variable u is 

u = L z L y L x u . ( 10 ) 
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In order to preserve second-order accuracy in space, in successive timestep the order of directional integration is permuted 
as follows: 


L z LyL x 


L x LyL z 


L x L z Lit 


L y L z L x 


LyL x L z 


L Z L X L„ 


( 11 ) 


For cosmological applications several important modifications must be introduced into the basic PPM scheme. The inclusion of 
gravitational forces and cosmic expansion changes the form of the usual hydrodynamical equations. Expansion of the universe 
enters the hydrodynamical equations in two different ways. First, all of the fluxes and the source terms are multiplied by the 
factor l/a(t). Second, two further terms appear: ( a/a)gvj and ( a/a)E in the equations for conservation of momentum and 
energy, respectively (in the following we will refer to these terms as expansion terms). This leads to important changes both 
in the Riemann solver, as we will show in the next section, and in the final integration. In this last step, forces and expansion 
terms are treated in different ways. First hydrodynamical quantities are updated without considering the expansion terms. 
These in fact, being proportional to the integrated quantity itself, cannot be included in the basic integration procedure. 
The time centered values of gravitational forces are calculated by finite differencing the potential field computed by linear 
extrapolation at t n+1 ' 2 . In order to obtain a proper estimate of the force felt by a fluid element, its averaged value over the 
entire zone is calculated. Finally, we compute the effect of expansion terms by using a semi-implicit scheme (see equations 
A8-A10 and A17). 


3.2 Characteristic Equations 


The solution of the Riemann problem, which is the central feature of the Godunov approach, requires knowledge of the 
characteristic form of hydrodynamic equations. Because of our choice of coordinates comoving with the expanding universe, 
the finite difference expressions which we use in the Riemann solver differ slightly from those discussed in Woodward & 
Colella (1984). In this section we write the characteristic equations for the system (1)—(3), while in Appendix A2 we present 
the scheme used in the code for computing hydrodynamic fluxes. 

For the time integration we use a directional splitting procedure. Indicating by x the relevant space direction during the 
one-dimensional sweep, by v the velocity in the direction of the one-dimensional sweep and by u the velocity orthogonal to 
v, equations (l)-(3) can be rewritten in the following non-conservative form: 


dU „ 

w +A w +c 


where 


0 , 


u 


f g \ 

( v 

0 

0 

e \ 

( ° \ 

V 1 , A = - 

0 

V 

0 

V* , c = 

va/a-g 

u I a 

0 

0 

V 

0 

ua/a I 

\pj 

Vo 

IP 

0 

V / 

V 2 pa/a ) 


( 12 ) 


(13) 


In order to find the characteristic form of this system, we first solve the corresponding eigenvalue equation: 


det (A — A T) = 0 


(14) 


obtaining as eigenvalues: 


Ac = - 

a 


A + — 


V + c 


(15) 


where Ao is a double solution and c is the sound speed. The eigenvalues are used to compute the domain of dependence of a 
given fluid element. Note that in expanding coordinates the domain must be rescaled by the factor a(t) with respect to that 
for a non-expanding system. 

The corresponding left eigenvectors R are: 


loi = (0,0,1,0) 

(16) 

1 02 = (1,0,0,-1/c 2 ) 

(17) 

!- = (0,-f,0,i) 

(18) 

4 = (o,f,o,i). 

(19) 

The characteristic form of hydrodynamical equations is then: 


, au , .au , ^ „ 

R —7— —b R^gJ dt — 0 , 

(20) 

, au , , au , ^ 

R ~dt —4C dt — 0 , 

(21) 
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that is: 

IferfU + lfcC dt = 0, 

For our fluid, the characteristic equations are: 


alon g ^ = A ,. 


c 2 dg — dp — 2—p dt = 0 
a 


du - u dt = 0 

a 

along dx/dt = v/a, and 
dp ± gc dv + 


_ a < 

(a V 

2— p dz qc 

-v- g) 

a ' 

va ) \ 


dt = 0 


( 22 ) 


(23) 

(24) 


(25) 


along dx/dt = (v ± c)/a. 

Characteristic equations are used both for calculating the domain of dependence of a zone edge and in the construction 
of the Riemann solver where we take account, to second-order, of the gravity and expansion (see Appendix A2). The solution 
of the Riemann problem gives the required time-averaged values of the hydrodynamical variables which are used in the final 
integration step. 


3.3 Correction to the Energy 


In order for overall energy to be accurately conserved it is important to use the total energy density E directly as a dependent 
variable (see equation 3) rather than calculate it as the sum of the kinetic and internal energy densities (^f?v 2 and e, 
respectively). The value of e is however required in order calculate the thermodynamic quantities such as pressure and 
temperature and the original PPM code calculates e as the difference between E and igv 2 . However, there is a problem 
with this for cosmological simulations where highly supersonic flows are often present. In these situations, the ratio between 
the kinetic and internal energies can reach values up to 10 8 in which case the errors in calculating e by the standard PPM 
procedure become much larger than the quantity itself. 

In our code, we solved this problem by calculating both the total and internal energy per unit comoving volume. In 
comoving coordinates the equation for the internal energy density e = p/(y — 1) is the following: 


de 

dt 


1 d a 

-~-a—\ ev A ~ 2—e - 
a oxk a 


_j_ dv k 

aP dxk 


(26) 


Hence, at each timestep, both the total and internal energy equations are solved and the results combined according to the 
local conditions in the fluid. The value obtained from equation (26) has to be only used when the internal energy cannot 
be calculated correctly as the difference between total and kinetic energy. This is usually the case in expanding or weakly 
compressing (comoving) regions, that is in regions in which one of the following conditions holds: 

V-v > 0 , i|Vp| < rji . (27) 

P 


If neither of these conditions is verified, the internal energy is computed according to the following criterion: 


e = 



if E — /E > r )2 
otherwise , 


(28) 


where p\ and rj 2 are suitable parameters smaller than unity, fixed on the bases of numerical experimentation. Various tests 
have suggested the use of the following values: 771 = 0.3 and p 2 = 0.005. 


3.4 N—body Code 

The dynamics of the dark component is described using a PM N-body code (Hockney & Eastwood 1981). In order to solve 
equations (5)-(6) the standard second-order leapfrog method has been replaced by a second-order two-step Lax-Wendroff 
method (see Appendix Al). The leapfrog scheme requires a constant timestep, while hydrodynamics makes it desirable to 
have a variable time interval. Density and forces are calculated by using the cloud-in-cell interpolation scheme. 


3.5 Timestep 

In determining the new timestep at the end of each loop we impose several constraints. First of all, we require that the 
maximum variation of the expansion factor a during a timestep is less than 2%. We impose also that dark matter particles 
move no more than half a cell in a single step and that the baryonic dynamics satisfies the Courant condition: 
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At < min 


C c a(t) Ax 

[c + max(|u*|, \v v \, |v*|)] 


( 29 ) 


where C c is the Courant number and the minimum is computed over all cells. We set C c = 0.5. Finally, we require that in a 
single step the timestep increases by no more than 20% from its previous value. 

Our tests indicate that the expansion condition is important only in the first part of the simulation, when the amplitude 
of the fluctuations is still small. At later times the Courant condition becomes the dominant one. 


4 TESTS OF THE CODE 


Our code has been subjected to a series of numerical tests, ranging from the purely hydrodynamical to cosmological ones, 
in one and more dimensions. In each case conservation of energy has been checked. For the cosmological tests this check has 
been performed by using the cosmic energy equation which for an ideal gas in an expanding frame is (Peebles 1980): 

^-(E + W) + ~{2E + W) = 0, (30) 

at a 

where E is the total (internal plus kinetic) comoving energy and W is the gravitational energy of the gas. 

Equation (30) can be rewritten as: 


— [a(E + W)\+aE = 0 

and integrated with respect to the expansion factor giving: 


(31) 


d 2 E(a 2 ) — aiE(ai) + / E da = — — aiW(ai)] , 


(32) 


where ai and (i 2 are the values of the expansion factor at two different times. 
We checked energy conservation computing the quantity: 


R = 


a 2 -E(a 2 ) — ai-E(ai) + f° 2 E da 
-[a 2 W(a 2 ) -aiW(ai)] 


(33) 


For perfect energy conservation R should remain as unity. The results obtained in our tests are presented in the following 
subsections. 


4.1 One-Dimensional Shock Tube Test 

In Fig. 1 we present the results of a standard shock tube test. In this case we solve the hydrodynamical equations (l)-(3) 
without the source terms on the right hand side and with a = 1. Initially the gas is at rest (v = 0) and density and pressure 
show a central jump from g l = 1, p t = 1 to g r = 0.1, p r = 0.1. The calculation is performed using a grid of 128 zones, the 
initial time is ti = 1 and the final time, shown in the figure, is to = 20. 

The numerical results (open circles) are in very good agreement with the analytical ones (solid lines). In particular, shocks 
are typically resolved within one cell and contact discontinuities are represented in four or five zones. In order to reduce the 
diffusion of contacts we have used, for the density field, the discontinuity detection algorithm of Colella & Woodward (1984). 
We have not introduced artificial diffusion and/or numerical flattening since numerical ripples around steep gradients are 
either absent or very small and do not affect the solution during the whole evolution. 


4.2 Adiabatic Expansion 

After having tested the hydrodynamic part of the code alone, we have introduced expansion and gravity and have concentrated 
on tests describing physical situations that are expected in cosmological simulations. In all of the following tests we have 
assumed an Einstein-de Sitter pure baryonic universe with vanishing cosmological constant and we fixed the value Ho = 50 km 
s' 1 Mpc^ 1 for the Hubble constant. The comoving size of the computational box is L = 64 /i -1 Mpc. For this model our basic 
units of normalization are: 

. t 3H % , 1 ^ 

00 87rG to (67rG^) 1 / 2 ’ (3 ) 

where G is the gravitational constant. The expansion factor depends on the cosmological time as 

a(t) oc t 2 ^ 3 . (35) 
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As a first test we have simulated the adiabatic expansion of an unperturbed universe, starting at time ti = 0.01 from a 
homogenous distribution of matter with velocity Vi = 100 km s -1 and temperature T — 4.6 K. In this situation equations 
(2)-(3) are easily integrated giving: 

v = Vi(a,i/a) , T = Ti(a,i/a ) 2 , (36) 


where a, v and T are the expansion factor, the velocity and the temperature respectively, at time t. At the final time to = 1 
numerical and analytical solutions are compared and the errors calculated as: 


AT 



Av 


\V-num V an 
Van 


(37) 


where the subscripts num and an indicate, respectively, the numerical and analytical values of the variables. We have indicated 
with Aa the maximum fractional variation of the expansion factor in a timestep. Using Aa = 0.1 we have obtained AT = 
6 x 10 -3 and Av = 3x 1(P 3 . The errors decrease to AT = 8 x 1(P 5 and Av = 4 x 1(P 5 for Aa = 0.01. The value of the 
parameter Aa used in the simulations is fixed by considering the balance between accuracy and computational time. 


4.3 Zel’dovich Pancake 


In this test we describe the formation of a one-dimensional pancake. The test is particularly important for cosmological 
studies as all of the related physics (hydrodynamics, expansion, gravity) is included. In addition, it represents the evolution 
of an initial sinusoidal wave and so can be seen as a single-mode analysis of the full three-dimensional problem. 

At the initial time, corresponding to a^ = 0.01, a sinusoidal velocity field has been imposed on the computational grid. 
The amplitude of the initial perturbation has been chosen such that the caustic forms at a c = 0.5. At m the density and 
temperature fields are uniform. We set g = 1 and T = 100 K. The final time corresponds to ao = 1. The same test has been 
repeated using eight different numbers of zones ranging from N c = 8 to N c = 1024. 

During the first part of the evolution, we are in the linear phase and the numerical results can be directly compared with 
the analytical solution which applies for a non-collisional fluid since the effect of the pressure is negligible in this phase. This 
comparison gives an estimate of the accuracy of the code. Moreover, using results obtained with different numbers of zones 
we can study the convergence of the numerical solution. The error is computed as: 


A g = 


1 

K 


E 


I e(Xj) - QejXj) | 

Qe{Xi) 


(38) 


where g e is the exact solution. We found that at a ~ 0.1, when shocks are not present, the error is below 1% already with 
only 16 zones. 

The non- linear phase of the evolution is the most interesting. The strong shocks and large gradients that form are in fact 
a severe test for any Eulerian numerical code. The results at ao = 1 obtained using grids of 256 and 1024 zones are compared 
in Fig. 2. The only important difference between the two solutions is the height of the central density peak, but this was 
expected as a consequence of the different resolutions in the two cases. All of the other features the solution are in very good 
agreement between the two simulations showing the convergence of the numerical solution. Shocks are resolved in one zone 
and there are no numerical postshock ripples. 

In Fig. 3 we present the results of the energy conservation test. We have calculated at the final time the conservation 
parameter R (defined in equation 33) and the total gravitational (IT), kinetic ( Ekin ) and thermal ( Eth ) energy of the fluid 
when a various number of cells N c are adopted. The energy is plotted in the code units. The error in the energy conservation 
is of about 5% for N c = 32 and decreases to 0.1% for N c = 1024. The behaviour of the three energy terms shows that at least 
32 zones are needed for starting to see the convergence of the solution. The poor resolution of the 8 and 16 zone cases, in fact, 
makes it impossible to describe correctly the narrow central peak, which remains too diffused. The shallower central potential 
wells induce a smaller acceleration in the infalling matter. The compression is low and the transformation from kinetic to 
internal energy is inefficient and so produces a pressure which is too low to support the infall of the matter. However the 
collapse of the peak is prevented by the artificial effect of the coarseness of the grid. 


4.4 Multiple Pancake Formation. 

The last one-dimensional test which we have performed consists of simulating the evolution of a random initial density field. 
In this test different wavelengths evolve together; hence we can verify the ability of the code for describing shock interactions 
and the merging of small structures. 

The initial conditions are set by perturbing the density field according to a Gaussian random distribution characterized 
by a power-law spectrum of the general form 

P{k) = P 0 k n exp(-k 2 R 2 f ) , (39) 
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where Po is the normalization constant and n is the spectral index. The short wavelength cut-off at the scale Rf = 2 grid- 
points ensures that the results are not affected by the sampling of modes whose size is close to that of the resolution of the 
simulation. 

The normalization constant is fixed requiring that the one-dimensional variance a 2 at the final time to (corresponding 
to the expansion factor a o = 1) for linearly evolved perturbations is equal to unity on the scale R * = 10^ 1 L so that: 


a 2 (R t ) = 


P{k)W{kR„) dk = 1 


where the function W(kR *) is a one-dimensional top-hat filter: 


(40) 


W(kR,) 


sin(A:P„) 

kRt, 


(41) 


The initial time ti is chosen by imposing that, at ti, the largest density fluctuation is equal to unity, so that the system can 
be considered still evolving linearly. 

As an example, we show in Fig. 4 the results at the final time to for a simulation with spectral index n = 1 when N c = 256 
is used, i.e. with the same number of cells used in the simulation of the single pancake. Shocked regions, characterized by 
the presence of typical double peaked maxima in the temperature field and by the smoothing of steep negative gradients in 
the velocity field, are at temperatures ranging from 10 4 to 10 6 K. Underdense regions are instead at much lower temperature. 
Collapsing regions surrounding density peaks are, at the beginning, at low temperature, but heat up as soon as propagating 
shocks reach them. The errors in energy conservation are around 3%, slightly higher than the value obtained for the Zel’dovich 
pancake test. 

The different clustering properties of baryonic and dark matter and the multiple pancake formation have been studied in 
much more detail in Gheller, Moscardini & Pantano (1996) by using numerical simulations based on a one-dimensional version 
of the PPM+PM code presented here. In that paper power-law initial power spectra with — 1 < n < 3 have been simulated 
using 8192 zones and 16384 particles. The results show that at large scales the final distribution of the two components is 
very similar while at small scales the dark matter presents a lumpiness which is not found in the baryonic matter. 


4.5 Cold Dark Matter Universe 


As first checks of our complete fully three-dimensional code, we repeated the tests of the one-dimensional shock tube and of the 
Zel’dovich pancake formation. Using a 64 3 grid, the results (not shown here) are consistent with those obtained with a similar 
resolution in one-dimension. Finally we have simulated the evolution of a purely baryonic Einstein-de Sitter Universe with 
an initial cold dark matter (CDM) power spectrum. As in the previous tests, we fixed the Hubble constant to be Ho = 50 km 
s' 1 Mpc -1 and the comoving size of the computational box to be L = 64 h _1 Mpc. In this way our numerical results can be 
compared with those obtained by Kang et al. (1994), even if a different random realization of the initial conditions has been 
used. Perturbations are initially Gaussian distributed and the primordial density power spectrum can be written in the form: 

P(k) = P 0 k n T 2 (k) , (42) 

where T(k) is the CDM transfer function given by the Davis et al. (1985) fitting formula: 

T(k) = (1 + 6.8fc + 72.0fc 3/2 + 16.0fe 2 ) -1 . (43) 


In the standard CDM model, the primordial spectral index has the Zel’dovich value n— 1. 

The normalization constant Po is fixed such that the linear mass variance at the present time is unity in a sharp-edged 
sphere of radius Rg = 8 hr 1 Mpc: 


cr 2 (P 8 ) = 1^/ k n+2 T 2 (k)W§. H (kR 8 ) dk= 1 . (44) 

In the above equation, WrH{kR) = ( 3/kR)ji(kR ) is a top-hat window function and ji denotes the Ith-order spherical Bessel 
function. The velocity field at the initial time ti has been calculated in the linear perturbation approximation (Peebles 1980) 
as: 


2 V$ 

3 H(ti)a{ti) 


(45) 


We are interested in checking the convergence properties of our code, so we evolved the same initial conditions on three 
different grids with 32 s , 64 3 and 128 3 cells, respectively. The initial conditions are set up according to the previous prescriptions 
on the coarsest grid. The initial time, fixed by the condition that the maximum density fluctuation is unity, corresponds to a 
redshift z ~ 20. We run our simulations on a IBM RISC 6000/590. 

We have first analyzed some integral properties of the results at the final time. We have calculated the mean temperature 
(T), the average mean temperature {T) e = (Tq)/{q), the average temperature weighted by the density squared (T ) e 2 = 
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[tp] 

Table 1. The simulations. Column 1: the number of zones in each dimension N c ; Column 2: the mean temperature ( T ). Column 3: the 
average mean temperature (T) B . Column 4: the average temperature weighted by the density squared (T) 2 . Column 5: the r.m.s. of the 
density fluctuation field a. 


N c 

(T) 

(T) e 

(T) g 2 

(7 

32 

1.48 

6.02 

19.49 

2.17 

64 

2.02 

10.61 

28.63 

2.48 

128 

2.02 

11.91 

32.95 

2.77 


(Tq 2 )/(q 2 ) and the variance of the density fluctuation field a 2 = ((q 2 )/(q) 2 — 1). As in Kang et al. (1994), the calculation has 
been performed for the data rebinned to a 16 3 grid. The results are presented in Table 1. 

The average temperature for the rebinned data seems to have converged already in the 64 3 simulation. However, the 
behaviour of the rms density fluctuation suggests that a higher resolution is needed for the convergence of the density field. 
The variation of the density weighted temperatures between the 64 3 and 128 3 grids depends essentially on variations of density 
distribution. As already noticed by Kang et al. (1994) for all the eulerian codes considered in that paper, the values of a 
appear to converge from below. 

I 11 Fig. 5 we show a cell by-cell comparison of density and temperature for the rebinned data of the three simulations. We 
always observe a good correlation between the values of the densities obtained from simulations with different resolutions; the 
spread of the points is reduced as we compare simulations with the highest resolution. The same holds for the high temperature 
cells. The distribution at lower temperatures is due to the formation of substructures as we increase the grid resolution. The 
large number of points which accumulate in the lower part of the panel represents cells that in the lower resolution simulation 
are cooling according to the expansion of the Universe while in the higher resolution run are site of some structure. Also this 
effect decreases when the resolution is higher. 

I 11 the left and right columns of Fig. 6 we show the volume and mass weighted density histograms, f(N) and f(M) 
respectively, for the 16 3 rebinned data. The tendency to form higher density peaks as we increase the resolution clearly 
appears in these histograms. As we go from 32 3 to 128 3 the amount of mass contained in the high density regions increases. 
However, in relation to the convergence of the numerical results, it is encouraging to see that the density distribution both in 
volume and in mass is very similar for the 64 3 and 128 3 runs. 

The corresponding histograms for the temperature are shown in Fig. 7. Comparing the rebinned data from different runs 
we can observe higher gas temperatures when the resolution of the simulation is increased. We notice both an increase in 
volume fraction and in mass fraction of the gas at high temperature. The f(N) peak at low temperature is almost absent in 
the 128 3 simulation, but this is only due to the averaging procedure in constructing the rebinned data. Small scale structures 
are present in the high resolution run which could not appear in the simulations with a smaller number of cells and they have 
the effect of raising the average temperature in the rebinned data. 

The same effect is evident if we look at Fig. 8 which presents the same histograms for the original data of the 128 3 
simulation: a large fraction of the volume is still at low temperature. However, the mass fraction distribution confirms that 
about 70% of the mass is at temperatures higher than 10® K. 

In Fig. 9 we show two different slices of the simulation with 0.5 /i _1 Mpc thickness (1 cell). For each slice contour plots are 
presented for the density contrast 8 — q/{q) — 1 and the temperature. Structures appear well resolved and strong temperature 
gradients indicate shock positions. High density peaks are surrounded by high temperature regions which have been heated 
by the shock propagation. 

Finally in Fig. 10 the contour plots of the volume (left) and mass (right) fraction with given temperature and density 
are shown. We observe that most of gas is at high density and temperature. On the other hand most of the volume is at low 
density with either high or low temperatures depending on whether it has been crossed by shocks or not. 


5 CONCLUSIONS 

We have presented a new numerical code developed for studying the formation and evolution of cosmological structures in 
both baryonic and dark components. Collisional matter is treated as a fluid and the corresponding hydrodynamic equations 
are solved using the PPM scheme on a fixed Eulerian grid. We have described the changes to the basic method required by 
the cosmological applications. Particular care has been taken in including expansion and gravity in the Riemann solver and 
in the final integration step. This has required the calculation of the characteristic form of the hydrodynamic equations in 
expanding coordinates. A double formulation of the energy equation has allowed a proper treatment of the highly supersonic 
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flows common in cosmological simulations. The behaviour of the dark matter is described using a standard Particle Mesh 
N-body technique, modified to allow the use of a variable timestep, as desirable for hydrodynamics. The two components 
are coupled through the gravitational interaction and the gravitational field is calculated from the Poisson equation using an 
FFT procedure. 

We have presented a series of tests selected for their relevance in cosmological applications, paying attention both to the 
accuracy of the highest resolution results and to the convergence of the method when lower resolutions are used. 

The one-dimensional tests show that the code can reproduce properly the expected solutions, even when very low 
resolution is adopted. In particular we present the results of the shock tube test and of single and multiple pancake formation. 
The CDM test results can be compared with those presented by Ryu et al. (1994), Kang et al. (1994) and Gnedin (1995) 
showing a good agreement. All of these tests have demonstrated that our code can be considered a reliable and useful tool for 
cosmological studies. 
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FIGURE CAPTIONS 

Figure 1. The velocity v, the pressure p and the density g for the one-dimensional shock tube test at the final time to = 20. 
The numerical result obtained with a grid of 128 zones (open circles) is compared to the analytical solution (solid line). 
Figure 2. The velocity v, the temperature T and the density g for the one-dimensional Zel’dovich pancake test at the final 
time cio = 1. The results obtained with grids of 256 (open circles) and 1024 zones (solid line) are compared. 

Figure 3. The energy conservation parameter R (top left), the gravitational energy W (top right), the kinetic energy Ekin 
(down left) and the thermal energy Eth (down right) for the one-dimensional Zel’dovich pancake test are plotted as function 
of the number of cells N c . 

Figure 4. The velocity v, the temperature T and the density g for the one-dimensional multiple pancake test at the final 
time ao = 1. A primordial power-law power spectrum with spectral index n = 1 is simulated using a grid of 256 zones. 
Figure 5. A cell by-cell comparison of density g (upper row) and temperature T (lower row) for the 16 3 rebinned data at 
z = 0 for three simulations with the same initial conditions but different numbers of zones, as indicated by the subscripts. 
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Figure 6. The volume-weighted f(N) (left column) and mass-weighted f(M) (right column) histograms at z = 0 for the 
density g at z = 0 for the 16 3 rebinned data. The results obtained with different numbers of zones are shown: N c = 32 (upper 
row), N c — 64 (central row) and N c = 128 (lower row). 

Figure 7. The same as Fig. 6 but for the temperature T. 

Figure 8. The volume-weighted f(N) (upper panels) and mass-weighted f(M) histograms (lower panels) for density g (left 
column) and temperature T (right column) at z = 0 for the original unbinned data of the 128 3 simulation. 

Figure 9. Density (left column) and temperature (right column) contour plots of two different slices of 0.5 hr 1 Mpc thickness 
in the 128 3 run are shown. Density is normalized to the mean density while the temperature is in units of 10 6 K. The density 
contour levels correspond to 10^ -3 ^ 4 , while the temperature levels are 5^ 2 3 (solid lines) and 5 -2 * (dashed lines), where 
i = 1 , 2 ,... 

Figure 10. Contour plots of the volume fraction (left) and mass fraction (right) with given temperature and density. The 
results at z = 0 from the 128 3 simulation are shown. Contour levels correspond to 10 l,/4 , where * = 1,2,.... 


APPENDIX A: NUMERICAL SCHEME 


We here describe the main steps in the integration of the dynamical equations presented in the text. We used a modified 
PPM scheme to follow the evolution of baryonic matter and a PM N-body method for the dark matter. The same grid of 
constant mesh size A® is used for solving the hydrodynamical equations and for computing the total gravitational potential. 
The subscript / indicates zone centred values of grid quantities, while subscript (/ + 1/2) refers to values computed at the 
boundary between the j —th and (/ + l)-th zones. Finally, in the dynamical equations for the dark matter the subscript / 
refers to the /-th particle. Superscripts indicate the time level at which variables are calculated. The same time interval 
Af n+1 ' 2 = t n+1 —t n is used for integrating the dynamical equations of both components. The time-integration of the N-body 
dynamical equations is performed by a Lax-Wendroff scheme instead of the usual leapfrog scheme which would require a 
constant timestep. 

The main steps in the integration of the dynamical equations are the following: 

1 - At the time t n we have g^Mj, v bm,j and Pj f° r the baryonic matter. For the dark matter we have Hbmj and v 2> M,j- 
First we calculate the dark matter density g'oMj using a cloud in-cell mass assignment scheme. Then we compute the total 
peculiar gravitational potential </" on the grid by solving the Poisson equation using a standard FFT technique. We also 
extrapolate the gravitational potential to £ n+1 / 2 using the values at t n and t n_1 


ra+l/2 _ 


1 + 


A£»+l / 2 

2At 71 - 1 / 2 


,_iAt" +1/2 

2Af"- 1 


(Al) 


2 We start the hydrodynamical one-dimensional sweep. For a given direction we indicate by v the component of the velocity 
along the direction of integration and by u the component orthogonal to v. We then calculate the time averaged values of the 
hydrodynamical quantities at the zone interfaces Qj+ 1 / 2 , pj+ 1 / 2 , Vj+ 1/2 and Uj+ 1/2 using the PPM scheme, as described in 
Section 3.1, 3.2 and in Appendix A2. 

3 - We use the time averaged estimates calculated in point 2 to solve the one-dimensional hydrodynamical equations without 
the inclusion of the expansion terms (for simplicity in the following equations we have dropped the BM subscript): 


n+1 _ n Af n+1/2 ( Qj-l/2Vj-l/2 - Qj+l/2Vj+l/2\ 
Qh, j — Bj + 


jn+1/2 


A* 


n-\-1 
V H,j 


O n + 1 

4 !HA 


n n 
QjVj 


At n+1/2 f Qj-l/2v'*_ 1 / 2 +Pi-l/2 - Qi+l/2Vj + 1 /2-Pj+l/2\ | Qj +1 + 0 j 

4 7, 9 3 


a 


n+ 1/2 


Aa; 


- 


(A2) 

(A3) 


£"+/ = £"+ 


At n+ 1/ 2 

a n+i / 2 


l/2V 2 j-l/2 + ^rTPj-l/ 2 ) Vj-1/2 ~ (lgj + l/2«| + l/2 + 7^iPj + l/2) Vj+1/2 

Ax 


n+l v n+ 1 + n v n 

+~ - 3 —z - —9j( A4) 


n+1 

' X H,j 


n+i 
9 HA 


n n . 

63 u j 4 “ 


At" + 1/ 2 / gj-l/ 2 Vj-l/ 2 Uj-l /2 — gj+l/2Vj+l/2Uj+l/2 


n+1/2 


Ax 


(A5) 


The quantity gj is the time-centred value of the gravitational force averaged on the /-th cell calculated by expanding the 
density and acceleration fields and retaining all terms up to second-order in Aa;: 


9 3 


2a(t)Ax 


<t>3 +1 ~ 0j-i + J2^ j+1 ~ ^ ^- 1 ) ^ 


Qi 


(A6) 
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The monotonicity of the numerical method is preserved by computing the density slope Sq as: 


$Qj = 0minSign(Q j+1 - Qj-i) ifpj+i - Qj){pj - Qj-i) < o 

= 0 otherwise , (A7) 

where Qmin = min(|gj+i — Qj-i\,2\6j ~ Qj-i\C\Qj+i ~ Qj\ )■ Because of its smoothness, the gravitational potential does not 
require the same monotonicity constraint. 

Points 2 and 3 are repeated for the other two orthogonal directions completing the three-dimensional hydrodynamical 
calculation. 

4 - The hydrodynamical quantities are corrected by considering the expansion terms: 

(AS) 


71 + 1 71 + 1 

Qi = en,i 


71 + 1 


e: 


i 71+1 


2v n H P - Af n+1/2 (d/a) n+1/2 < 
2 + At"+ 1 / 2 (a/a)"+ 1 / 2 

E n H +f - At n+1/2 (d /a) n+1/2 E? 
1 + A t n + 1 / 2 (a/a) rl + 1 / 2 


(A9) 

(A10) 


5 - Position and velocity of the N-body particles are updated by the two-step Lax-Wendroff scheme: 
step 1 

Af n+1/2 


n+l/2 
* DM, j 


n+l/2 _ 

step 2 


71 ( • 71 71 i 71 \ 

— v Btf,j — p v DM,j + Sj ) 

a r+ 1/2 


2 a n 


^DM,] + v DM,j‘ 


2a" 


v «+! — v - 

v DM,j — v DM,j 


_ ( ^+1/2 n + l/2 n+l/2\ At" + 1/2 

1 y a V DM,j -T&j J a n + l/2 


71 + 1 

DM,j 


, Af n + 1 / 2 

n+l/2 AAT_ 


— x DM,j + V DM,j n+l/2 


(All) 

(A12) 

(A13) 

(A14) 


Here we indicate by gj the gravitational force on the j—th particle, computed by interpolating the values of the corresponding 
components at the eight neighbouring cells using the cloud-in-cell scheme. 

6 Energy is corrected to account for highly supersonic flows. We have introduced the internal energy equation. This is solved 
first including the flux term by the standard PPM integration procedure: 

A t n+1/2 f Pj-l/2Vj-l/2 ~ Pj+l/2Vj+l/2 


~71+1 /t , 

e H,j = e* + 


2 a(t n+1 / 2 ) 

Then this value is corrected as: 

1 At n+1/2 P- +1/2 


(7 — l)Ax 


n+i 

e H,j 


~n+l 

e HJ 


i +!/2 


{Vj- 1/2 - Vj+ 1 / 2 ) 


where pP 1 ^ = (pj - 1/2 PPj+ i/ 2)/2. Finally, we consider also the expansion term in the internal energy equation: 

n+1 = e H + j - At n+1 ' 2 (a/a) n+1 ' a e? 

6j 1 + At n+1 / 2 (a/a) n+1 P 

The internal energy is used if conditions (27) hold: 

1/2 

'dp 


E 

7=1,3 


dvj 

dxi 


> 0 


1 

Pj 


e (m 


< 771 j=l,3. 


If neither of these conditions are verified, internal energy is computed according to expression (28): 

Ej - \QjV 2 if Ej - \QjV 2 1 max(£j_i, Ej, E j+1 ) > 
e-j otherwise . 


(A15) 


(A16) 


(A17) 


(A18) 


(A19) 


Note that as we are dealing with errors that can be advected with the total energy calculation, we have to consider not just 
the j-th cell but also its neighbourhood (that in one dimension comprises the two zones j ± 1). 
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APPENDIX B: CALCULATING FLUXES 


Here we describe briefly the procedure used to calculate the time averaged values of the hydrodynamical quantities at the 
zone interfaces Qj+ 1 / 2 , Pj+ 1 / 2 , iij+ 1/2 and Vj+ 1/2 • We include to second-order both the gravity and the cosmic expansion, 
which enter the characteristic equations as described in Section 3.2. 

We indicate by w(x) the interpolated piecewise continuous function representing the generic hydrodynamical variable. 
The quantity Wj represents the mean value of w(x) over the ji—th zone. 

We define u>j+i/2,L and Wj+i/2,R, first approximations to the effective left and right states of the Riemann problem, as: 


V>j+1/2,L — fj+1/2 ,l( x j+1/2 — Xj+1/2,L,) , 
V>j+1/2, R = /j + l/2,Ji( — *j + l/2 + *j + l/2,J?) 

where 


/i+i/2f( Al ) ^ 


x j+l/2 


«’(£) d £, 


' x j + l/2~^ x 

Xj+l/ 2 ,l = x j+1 /2 - max(0, At n+1/2 Aj,+) 
and 


/j+1/2 ,h(A:e) - — 


x j+l/2 


x i + 1/2 


+A.x 


«’(£) d£ 


Xj+1/2 ,r = Xj+l/2 +max(0,-Ai" +1/2 Aj+i,_) 


-\j,+ — 


A /+i,- = 


u j +1 


U'+i 


(Bl) 

(B2) 

(B3) 

(B4) 

(B5) 


a"- a"- 

We correct our initial guess for the left and right states by solving the equations of gas dynamics in characteristic form 
(see Section 3.2). First we calculate the mean value of the hydrodynamic variables in the domains defined by each of the 
characteristic lines coming from the left and the right of the zone boundary j + 1/2: 


w j+l/2 ,L — fj+l/ 2 ,L{Xj+i /2 ~ *J+1/2,l) j 


(B6) 

(B7) 


w j+ 1 / 2 ,R ~ fj+ 1 / 2 ,r(— Xj+ 1/2 + Xj+i/ 2 ,r) , 
where k = +, 0, — and 

Xj+1/2, l — Xj+ 1/2 xj+ 1 / 2 ,r ~ Xj+i /2 AfAj-j-ufc . (B8) 

Then we correct the initial guess considering only that information which can reach the zone edge during the current timestep: 
Vj+l/ 2 ,S = Pj+1/2,S + Cj+ 1 / 2 , s{Pj+l/2,S + Pj+l/2,s) > ( B 9) 

rj+ 1 / 2 ,S = Vj+ 1 / 2 ,s + C j+l/ 2 ,s{h++l/ 2 ,s ~ Pj+ 1 / 2 ,s) > (BIO) 

( 1 


£b+i /2 ,s 


Pj+1/2,S 


- E s 

fc=+,0,- 


j + 1/2,5 


(Bll) 


Here Cj +1 / 2 ,s = 7 Pj+i/ 2 f?j+i /2 and S = L, R. We also have 


V j+1/2,L — Vj+1/2,L, 
Vj + l/ 2 , R = Vj+i/ 2 ,R, 

Otherwise 


Pj+1/2 ,l — 0 if A j,k < 0 , 
Pj+1/2 ,r ~ 0 if A j+i,fc < 0 • 


P: 


,± 

0 + 1 / 2 , s 


1 


C 


j+1/2,s 


(Vj+1/2,S - vf + i/ 2 ,s) ± 


Pj+1/2, S Pj+1/2, S 


C 


± At 2- 


a Pj+ 1 / 2 ,s 


j+1/2, s 


aC. 


j+1/2, s 


Pj+1/2,S - 


Pj+1/2, S - Pj+1/2, S 


1 


r 2 

J + 1/2,s 


Pj+1/2,S ^'+ 1 / 2,3 


+ q « P j+1/2,S 
a ^j+1/2, S/ 


Vj+l/2,S — Vj + i/2,s 



(B12) 

(B13) 

, 0 ± ± \ 

± 0 V j+l/2,S T Pj + 1/2,S 1 

(B14) 


(B15) 


(B16) 


The time averaged estimates of the solution at Xj+ 1/2 are determined by calculating the non-linear interaction of the left and 
right states calculated in equations (A33)-(A35). This interaction is described by the solution to the Riemann problem with 
these left and right states as initial conditions. The Riemann problem is solved in a standard way, as described by Van Leer 
(1979) and Gottlieb & Groth (1988). 
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